投影数据预处理与重建

说明

数据

​ ParallelBeam_ProjectionDara.raw ：样品投影数据。采用平行束CT系统对样品进行CT扫描，探测器的通道数为512，在360度范围内等间隔扫描，扫描角度间隔为1度，共扫描了360个角度的投影数据，数据大小512×360。

​ ParallelBeam_DetectorFlatData.raw：探测器平场数据（亮场数据）。该平场数据仅在一个特定角度采集，采集时不放扫描样品，数据大小512×1；两数据均采用“unsigned short”格式保存。另外已知该探测器存在少量坏点，探测器的暗场影响忽略不计。

要求

  1. 请对该系统采集的投影数据作预处理，进行平场校正、坏点补偿，以及对数运算操作，然后利用滤波反投影重建算法进行断层重建。

  2. 请撰写实验报告，详细说明数据处理流程，展现各个处理步骤的结果，同时展示和分析各个步骤对重建结果的影响，并附上代码。实验报告要求格式清晰，讨论分析清楚充分，写明姓名学号，提交作业时，文件命名“5_PreProcess_学号_姓名.doc” 或“5_PreProcess_学号_姓名.pdf”。

方法

1. 数据预处理

使用fopen和fread读取裸数据。

原始投影数据 
% 读取样品投影数据
f1 = fopen('ParallelBeam_ProjectionDara.raw','rb');
I_sample = fread(f1,[512,360],'ushort');
fclose(f1);

% 读取平场数据
f1 = fopen('ParallelBeam_DetectorFlatData.raw','rb');
I_bright = fread(f1,[512,1],'ushort');
I_bright = repmat(I_bright,1,360);
fclose(f1);

2. 平场矫正

对投影数据进行矫正

Isample_dark=Isample−IdarkI_{sample\_dark} = I_{sample}-I_{dark}

Ibright_dark=Ibright−Idark I_{bright\_dark} = I_{bright}-I_{dark}

Isample_sample_bright=Ibright_mean∗Isample./Ibright_darkI_{sample\_sample\_bright} = I_{bright\_mean}*I_{sample}./I_{bright\_dark}

其中IsampleI_{sample}为扫描的投影图；IbrightI_{bright}为扫描的平场数据；IdarkI_{dark}为暗场数据。在本文中忽略暗场数据，设为0。Ibright_meanI_{bright\_mean}为$I_{break_dark}的均值。

平场矫正 
% 暗场矫正 直流偏置矫正
I_dark = 0;
I_sample_dark = I_sample- I_dark;
I_bright_dark = I_bright- I_dark;

% 平场矫正 通道一致性矫正
I_bright_mean = mean(mean(I_bright_dark));
I_sample_dark_bright = I_bright_mean*I_sample_dark./I_bright_dark;

3. 坏点补偿

在原始图像中可以看到投影数据中的坏点，通过对数据进行筛选，确定坏点为探测器编号100和220。使用线性插值的方法补齐数据。

坏点补偿 
I_sample_dark_bright_deBL = I_sample_dark_bright;
I_sample_dark_bright_deBL(100,:) = BedLine(I_sample_dark_bright,100);
I_sample_dark_bright_deBL(220,:) = BedLine(I_sample_dark_bright,220);

function result = BedLine(Img,number)
% 用于计算坏点数值
result = 0.5*(Img(number-1,:) + Img(number+1,:));
end

4. 对数运算

X光射线穿过被摄物体到探测器，符合衰减函数

I=I0e−udI= I_0e^{-ud}

转换得到图像的投影信息

prjpre=−ln(I/I0)prj_{pre} = -ln(I/I_0)

对数计算 
I_sample_dark_bright_deBL_deSA = DeSignalAttenuation(I_sample_dark_bright_deBL);

function data_log = DeSignalAttenuation(data)
% 对接收信号做衰减逆处理
I0 = max(max(data));
data_log = -log(data/I0);
end

5. 滤波反投影重建

使用iradon函数计算，之前几节已经讨论过重算算法，这里只展示结果。

重建图像 
Img = iradon(I_sample_dark_bright_deBL_deSA,1:360,'linear','Ram-Lak',1,512);

结果分析

步骤比较
重建结果比较

通过比较以上四幅不同步骤下的图像，可以看到：

  1. 平场矫正使得投影数据更加均匀，消除探测器噪声重建图像较好。

    平场矫正切线

    取投影图180度的位置绘制切线图，第一张为原始数据，第二张为做了平场矫正数据。可以看到平场数据信号较为平缓，消除了原始信号存在的高频小幅噪声。在坏点方面，将原始图像的杂乱信号标准化为探测器的基础信号，方便下一步处理。

  2. 坏点补偿使得信号更为连续，重建图像对比度高较为清晰。

    坏点补偿切线

    可以看到之前缺失的坏点信息已经补上了。

  3. 对数运算主要操作，用于将接收到的信号转化为图像重建当中具有物理含义的信号，表示为穿过体素的密度。对应关系变化，背景信号处理为零。

对数运算切线